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Abstract 

A self-adaptive moving mesh method is proposed for the numerical simulations 
of the Camassa-Holm equation. It is an integrable scheme in the sense that it 
possesses the exact A^-soliton solution. It is named a self-adaptive moving mesh 
method, because the non-uniform mesh is driven and adapted automatically by 
the solution. Once the non-uniform mesh is evolved, the solution is determined 
by solving a tridiagonal linear system. Due to these two superior features of the 
method, several test problems give very satisfactory results even if by using a 
small number of grid points. 
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1. Introduction 

Since its discovery [HI], the Camassa-Holm (CH) equation 

wt + 2kwx - wtxx + "iwwx = 2wxWxx + wwxxx- (1) 

has attracted considerable interest because it describes unidirectional propagation 
of shallow water waves on a flat bottom. It also appeared in a mathematical search 
of recursion operators connected with the integrable partial differential equations 
12] ■ By virtue of asymptotic procedures, the CH equation was reconfirmed as 
a valid approximation to the governing equation for shallow water waves Isl, 0]. 
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The CH equation also arises as a model for water waves moving over an under- 
lying shear flow [7], in the study of a certain non-Newtonian fluids [8], and as a 
model for nonlinear waves in cylindrical hyperelastic rods [9] . The CH equation 
is completely integrable (see [IJ for the Lax pair formulation and [fiol . [Till for the 
inverse scattering transform), and it has various exact solutions such as solitons, 
peakons, and cuspons. When k = 0, the CH equation admits peakon solutions 



which are represented by piecewise functions Illj,l3|,|4|]. When k 7^ 0, cusped soli- 
ton (cuspon) solutions, as well as smooth soliton solutions, were found by several 

authors. [|i3,[i3L[ii,[ii,H[i3,[i8L[i9]. 

Several numerical schemes have been proposed for the CH equation in the 
literature. These include a pseudospectral method 1I20I1 . finite difference schemes 
112 iL l22n. a finite volume method [23], finite element methods ll24l l25l. 12611. multi- 
symplectic methods [27], and a particle method in terms of characteristics based 
on the multi-peakon solution [28, .29.. i30 . 31. 32ll. We comment that the schemes in 
rt2lL I22I1 and in [1271] can handle peakon- antipeakon interactions. However, it still 
remains a challenging problem for the numerical integration of the CH equation 
due to the singularities of cuspon and peakon solutions. 

In the present paper, we will study an integrable difference scheme for the CH 
equation fll) based on an integrable semi-discrete CH equation proposed by the 
authors 1 3 31]. The scheme consists of an algebraic equation for the solution and 
the non-uniform mesh for a fixed time, and a time evolution equation for the mesh. 
Since the mesh is automatically driven and adapted by the solution, we name it a 
self-adaptive moving mesh method hereafter. 

As a matter of fact, Harten and Hyman has proposed a self-adjusting grid 
method for one-dimensional hyperbolic problems [34]. Since then, there has been 
significant progress in developing adaptive mesh methods for PDEs [|35l . 



38Ll39|,|40|]. These methods have been successfully applied to a variety of physical 



and engineering problems with singular or nearly singular solutions developed in 
fairly localized regions, such as shock waves, boundary layers, detonation waves, 
etc. Recently, an adaptive unwinding method was proposed for the CH equation 
II23I1 . The method is high resolution and stable. However, in order to achieve a 
good accuracy, a large number of grid points (= 4096) has to be used. In addi- 
tion, the designed method is only suitable for the single peakon propagation and 
peakon-peakon interactions, not for the peakon-antipeakon interaction. As shown 
subsequently, the self-adaptive moving mesh method gives accurate results by us- 
ing a small number of grid points 100) for some challenging test problems. 

The remainder of this paper is organized as follows. In Section 2, we present 
the self-adaptive moving mesh method and show it is consistent with the CH equa- 



2 



tion as the mesh size approaches to zero. Two time advancing methods in imple- 
menting the self-adaptive moving mesh method are presented in Section 3. In 
Section 4, several numerical experiments, including the propagations of "peakon" 
and "cuspon" solutions, cuspon-cuspon and soliton-cuspon collisions, are shown. 
The concluding remarks are addressed in Section 5. 



2. A self-adaptive moving mesh method for the Camassa-Holm equation 



It is shown in 113 311 that the CH equation can be derived from the bilinear equations 
of a deformation of the modified KP hierarchy 



- Qa/^x - 1 

2cff={D, + 2c)g-h, 
^ -2ff = {DtD^ + 2cA 

through the hodograph transformation 



f-f = 9h, 

2)g-h, 



(2) 



X = 2cx + log f 

r = t, 



(3) 



and the dependent variable transformation 



w 



log 



h 



f{x)g{y) 



y=x 



Here c = 1/k, and Dt are Hirota's D-operator defined as 

d__d_ 

dx dy ^ 

It is proved in that the bilinear equations ^ admit a determinant solution 
f = "Tq^ 9 = T-i, h = Ti, where r„ is a Casorati-type determinant of any size. By 
discretizing the x-direction with an uniform mesh size a, the following bilinear 
equations 

-Dt — 1 ) fk+i ■ fk = 9k+ihk + 9khk+i , 
a J 

'^acfk+ifk = (1 + 0'C)gk+ihk - (1 - ac)gkhk+i , 

-2afk+ifk = ((1 + ac)Dt - a)gk+i ■ hk - ((1 - ac)A + a)gk ■ hk+i , 



(4) 



3 



admits Casorati-type determinant solution with discrete index which is presented 
afterwards. Starting from Eq.©, a semi-discrete CH equation 



Wk+l -Wk Wk- Wk-1 



dSk 
dt 



1-^ 



5, 



fc-i 



Wk+l +Wk , 4 
Ok \ 



1 



SI 



+4-1 

{Wk+l - 



Wk + Wk-1 , 4-1 



1 - 



c 1 — a^c^ 



(5) 



0, 



Wk 



was proposed (see the details in 113311 ). Here the solution t) is approximated 

by Wk{t) at the grid points Xk (k = 1, ■ ■ ■ , N). The mesh 4 = Xk+i — Xk is a dis- 
crete analogue of the hodograph transformation from the x-domain with uniform 
mesh size a to X-domain. As is seen, it is non-uniform and time-dependent. 
The semi-discrete CH equation (|5]) can be rewritten as 



A'^Wk = —M 5kMwk + 



d&k 
dt 



c6k — a? 



(6) 



SkAwk 



by introducing a forward difference operator and an average operator A and M 

Fk+i — Fk , , _ Fk + Fk+i 



4 "2 
In the present paper, Eq. ^ or Eq. © is used as a numerical scheme for the 
CH equation ([T]). It is shown to be integrable in in the sense that it possesses 
A^-soliton solution which, in the continuous limit, approaches A^-soliton solution 
of the CH equation. The A^-soliton solution is of the form 



Wk 



log 



9k_ 
hk 



(7) 



with 



fk = ro{k) 



9k = ri{k), hk = T_i{k) 



(n+1) 



4"^ 



(n+Af-l) 



{n+N-l) 
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where 

i^t^ = a^fe - c)"(l - api)-^e^^ + a,,2(-p. - c)"(l + ap,)-'=e^' , 

= 1 + 4io , Vi = —i + Vio ■ 

Pi - c Pi + C 

Next, let us show that in the continuous limit, a — )■ (5,^ — )■ 0), the proposed 
scheme is consistent with the CH equation. To this end, the equation ^ is rewrit- 
ten as 



( -2 SkMwk , 6k-iMwk-i , 1 
[Awk - Awk-i) + - — — h 



h + h-i h + h~i h + h-i c(l - a'^c^ 

Aa'^c 1 



1 - a2c2 45fc_i 



=(1-^1 (^fc+i - U)k) ■ 



4 

By taking logarithmic derivative of the first equation, we get 



^* 1 :r^A — (^^^ ~ ^Wk-i) - ^— T — 

{Ok + Ok-l Ok + Ok^l Ok + Ok^l 

2 hMwk Sk-iMwk-i 

[Awk - Awk-i) 



5k + 6k-i 6k + Sk-i Sk + Sk-i c(l - a^c^) 

^k-l 



dt5k = ( 1 - ^ ) {wk+i - Wk) 



Thus, we have 



, 2 ^ 5kMwk 6k-iMwk^ 



4 + 4-1 Sk + 4-1 4 + 4-1 

2 6kMwk 6k-iMwk-i 

[Awk - Awk-i) 



4 + 4-1 4 + 4-1 4 + 4-1 c(l - a^c^) 

-{!-$] Awk-(l-^-^] Awk-i 
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The dependent variable w is a function of k and t, and we regard them as a 
function of X and T, where X is the space coordinate of the k-th lattice point and 
T is the time, defined by 

fc-i 

X = Xo + ^5,, T = t. 

j=0 

Then in the continuous limit, a — )■ (5^ — > 0), we have 

Awk wx , Awk-i wx , Mwk w , Awk^i -)■ wx , 

and ^ 

— — - — {Awk - Awk-i) wxx- 

Ok + Ok-l 

Further, from 



k-l 



dX dX, J^d5, dX, 5|\ 



at 

we have 



j=0 



ax 



a* = St + ^dx dT + wdx 
at 



where the origin of space coordinate Xq is taken so that 



dXo 
dt 



cancels wq. Then 



the above semi-discrete CH equation converges to the CH equation 

{dr + wdx){wxx - w) 



-2wx 



Wxx -w 

c 



I.e. 



((9t + wdx){wxx - w) = -2wx ( wxx - w 



(8) 



Setting c = we obtain the the CH equation ©. 

Note that, in our previous paper 113 311 . we put c = 1/ which gives an alterna- 
tive form of the CH equation 



wt + 2k Wx — wtxx + "iwwx = 2wxWxx + wwxxx- 



(9) 



It is shown that they are equivalent under the scaling transformation w — )■ kw, 
T — 7- T/k. In the present paper, for the convenience in comparing of our results 
with other papers QHlMoHH set c = 
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3. Implementation of the self-adaptive moving mesh method 



In this Section, we will discuss how to implement the self-adaptive moving 
mesh method in actual computations. Generally, given an arbitrary initial condi- 
tion w{X, 0) = Wo{X), the initial non-uniform mesh Sk can be obtained by solv- 
ing the nonlinear algebraic equations by Newton's iteration method. However, for 
the propagation or interaction of solitons or cuspons, which are challenging prob- 
lems numerically, the initial condition Wk can be calculated by (|7]) from gk and 
hk by putting t = 0, which are obtainable from the corresponding determinant 
solutions. The initial non-uniform mesh 5° can also be calculated by B3ll 

^0 ^ o (1 + ac)gk+ihk - (1 - ac)gkhk+i ^^^^ 
^ (1 + ac)gk+ihk + (1 - ac)gkhk+i ' 

On the other hand, once the non-uniform mesh 5k is known, the solution Wk 
can be easily obtained by solving a tridiagonal linear system based on the first 
equation of the scheme. 

a;<_Y + hiw^^^ + Q<Y = du (11) 

where 



and 



d, 



1 - a2c2 \5l+^ 51+1 ) c(l - a2c2) 



In regard to the evolution of 5k, we propose two time advancing methods. 
The first is the modified forward Euler method, where we assume Wk remains 
unchanged in one time step. Integrating once, we have 

where = (2 + 5^)/ [2 — 5"^). The second is the classical 4th-order Runge- 
Kutta method, where Wk can be viewed as a function of 5k by solving the above 
tridiagonal linear system. Therefore, in one time step, we have to solve tridiagonal 
linear system four times. 
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In summary, the numerical computation in one time-step only involves a ODE 
solver for non-uniform mesh and a tridiagonal linear system solver. Hence, the 
computation cost is much less than other existing numerical methods. A Matlab 
code is made to perform all the computations. Iterative methods, for instance, 
the bi-conjugate gradient method bicg in Matlab are used to solve the tridiagonal 
system. 

For the sake of numerical experiments in the subsequent section, we list exact 
one- and two- soliton/cuspon and peakon solutions. 

(1). One soliton/cuspon solution: The r-functions for the one soliton/cuspon 
solution are 



g (xl± 



c-Pi 
c + pi 



h(xl± 



c-pi 



(13) 



with ^1 = pi{2x — vit — xio), f 1 = 2/(c^ — pf). This leads to a solution 

2plcvi 



w{x, t) 



(c^ +p\)±{c^ — p\) cosh ^1 



X = 2cx + loE 



T = t. 



(14) 



(15) 



where the positive case in Eg. ([14)) stands for the one smooth soliton solution when 
Pi < c, while the negative case in Eq. (fT4l) stands for the one-cuspon solution when 
Pi > c. Otherwise, the solution is singular. Thus Eq.(fT4l) for nonsingular cases 
can be expressed by 



w{x, t) 



2p\cvi 



(c^ + pf) + — p\ \ cosh,^! 
Similarly, for the semi-discrete case, we have 

k 



(16) 



oc 1 



c-pi 



c + pi 



1 + api 
1 — api 



hk ccl + 



c + pi 



c-pi 



1 + api 
1 — api 



with ^1 = pi(—vit — Xio), resulting in a solution of the form 

2plcvi 



Wk{t) 



(C2 + pI) + ^ 



/ l+gpi 



=-6 



l+gpi 
1-api 



=5i 



(17) 



(18) 
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in conjunction with a transform between an uniform mesh a and a non-uniform 
mesh 

^ _ ^ (1 + ac)gk+ihk - (1 - ac)gkhk+i 
^ (1 + ac)gk+ihk + (1 - ac)gkhk+i ' 
(2). Two soliton/cuspon solutions: The r-functions for the two soliton/cuspon 
solution are 



5' ocl + 


c 


-Pl 


+ 


C-P2 


+ 


c 


+ Pl 




C + P2 




/l oc 1 + 


c 


+ Pl 


+ 


C + P2 


e^^ + 




c 






C-P2 





(C-Pi)(c-P2) 



(C + Pi)(c + P2) 
(C + Pi)(c + P2) 



(C-Pi)(c-P2) 



Pi + P2 



Pi +P2 



with ^1 = pi{2x - vit - xio), ^2 = P2(2a; - - X20), Vi = 2/{c^ - pi), 
V2 = 2/ [c^ — p\). The parametric solution can be calculated through 



(log|)^, X = 2cx + log(|), T^t, 



(19) 



whose form is complicated and is omitted here. Note that the above expression 
includes the two-soliton solution (pi < c, p2 < c), the two-cuspon solution (pi > 
c, P2 > c), or the soliton-cuspon solution (pi < c, p2 > c). 
Similarly, for the semi-discrete case, we have 



oc 1 + 



c-pi 



+ 



c + pi 
c- Pi){c- P2) 



l-apij 



C-P2 



1 + ap2 
1 - ap2 



(C + Pi)(c + P2) 



Pi -P2 
Pi +P2 



C + P2 

1 + api 
1 — api J V 1 ~ ^P2 



^ ^ 1 + ap2 ^ 



/ife oc 1 + 



c + pi 



c-pi 



1 + api 
1 - api 



C + P2 



C-P2 



1 + ap2 
1 - ap2 



={2 



+ 



(C + Pi)(c + P2) 



(C-Pl)(c-P2) 



P1+P2J \l-api) \l-ap2J 



with ^1 = pi{—vit — xio), {2 = P2{—V2t — X20). The solution can be calculated 
through 



w{x, t) = I log 



9k_ 
h 



(20) 



with a transforai 



[1 + ac)gk+ih 



ac)gkh 



k+l 



[1 + ac)gk+ihk + (1 - ac)gkh 



(21) 



fc+i 



Again, the explicit form of the solution is complicated and is omitted here. 

(3). Peakon solutions: In the continuous CH equation, it is possible to con- 
str uct p eakon solutions from soliton solutions by taking the peakon limit 41, 

For the continuous case, we can express the 1 -soliton solution as 



w 



1 + p{h? + (1 — K^)cosh^i 



where K = \,vi = 2k^/(1— /t^), = ]9i/t(2x//t—(fi/K)t—a;io/K). Taking the 
peakon limit /t — )■ 0, piK l,vi = const., the solution {X{x, t),w{x, t)), where 
X{x, t) = 2x/k + log |, gives the 1-peakon solution [l43[|. In FiglH one can see 
that the 1 -soliton solution approaches to the 1-peakon solution as k approaches to 
0. 





Figure 1: 1-soliton solution for the CH equation: the left: pi = 0.5, c = 1; the right (close to the 
peakon limit): pi = 99, c = 100. 



We can also consider the peakon limit for the semi-discrete CH equation. For 
the semi-discrete case, we can express the 1-soliton solution as 



1 + pIk^ + 



l—p^K? 



y 1-api 



l+gpi 
1-api 



=6 



where k = ^, vi = 2k^/{1 — p\h?), = pik{—{vi/ K)t — xiq/k). The peakon 
limitfor the semi-discrete CH equation is again K — )■ 0,pi/« — 1, t>i = const. Tak- 
ing the peakon limit, the solution {Xk{t),Wk{t)), where Xk{t) = -^o+Z]j=o ^ji^)' 
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approaches to a solution which approaches to the peakon solution of the CH equa- 
tion as taking the continuous limit. In Figl2l one can see that the 1-soliton solution 
approaches to the 1 -peakon like solution as k approaches to 0. Taking the contin- 
uous limit, this solution approaches to the 1 -peakon solution of the CH equation. 




-20 -li -ID -5 0^ -20 -10 10 



Figure 2: 1-soliton solution for the semi-discrete CH equation: the left: pi = 0.5, c — l^a — 0.1; 
the right (close to the peakon limit): pi ~ 99, c = 100, a = 0.005. 



4. Numerical experiments 

In this section, we apply our scheme to several test problems. They include: 
1) propagation and interaction of nearly-peakon solutions; 2) propagation and in- 
teraction of cuspon solutions; 3) interactions of soliton-cuspon solutions; 4) non- 
exact initial value problems. 

4.1. Propagation and interaction of nearly-peakon solutions 



Example 1: One peakon propagation. It has been shown in [|42l 14311 that 
the analytic A^-soliton solution of the CH equation converges to the nonanalytic 
A^-peakon solution when k — )■ (c — )• oo). To show this, we choose one soliton 
solution with parameters c = 1000, p = 998.9995. Thus the speed of the soliton 
(fi/2) is 1.0. Its profile is plotted and is compared with one peakon solution 
u{x,t) = e~l^~*l in Fig. |3l These two solutions are indistinguishable from the 
graph. The error in Loo, where = max \wi — ui\, is calculated to be O(10~^), 
and the discrepancy for the first conserved quantity Ii = J udx is less than 0.7%. 
Therefore, this soliton solution can be viewed as an approximate peakon solution 
with amplitude 1.0. 

The propagation of the above designed approximate peakon solution is solved 
by the self-adaptive mesh scheme with two different time advancing methods: 
the modified forward Euler method (MFE) and the classical Runge-Kutta method 
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Figure 3: Comparison between one peakon solution and one-soliton solution with c = 200.0. 



(RK4). The length of the interval in the x-domain is chosen to be 0.02 and the 
number of grid is = 101. For the above parameters of one-peakon solution, 
the length of the computation domain turns out to be about 28.5. Figures |4](a)-(d) 
display the numerical solutions at t = 1.0,2.0,3.0,4.0, together with the self- 
adjusted mesh. It can be seen that the non-uniform mesh is dense around the 
crest. The most dense part of the non-uniform mesh moves along with the peakon 
point with the same speed. With the same grid points = 101, the relative 
errors in Loo -norm and the first conserved quantity h = J w dx are computed 

and compared in Table \T\ Here, Loo = max , where wi and wi represent 

the numerical and analytical solutions at the grid points Xi, respectively. Ei = 
\h — h\/\h \ indicates the relative error in Ji, where Ji stands for the counterpart 
of Ii by the numerical solution. Trapezoidal rule on the non-uniform mesh is 
employed for the evaluation of the integrals. 

Example 2: Two peakon interaction. For c = 1000, we initially choose two 
approximate peakon solutions moving with velocity vi/2 = 2.0, and V2/2 = 1.0, 
respectively. Their interaction is numerically solved by MFE and RK4, respec- 
tively, with a fixed grid number of = 101. Figure [5] displays the process of 
collision at different times. Table [21 presents the errors in Loo -norm and Ei. It 
could be seen that, in spite of a small number of grid points and a large time step, 
RK4 simulates the collision of two approximate peakons with good accuracy. 

In regard to the propagation and interaction of approximate peakon solutions, 
we summarize as follows: 
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Figure 4: Numerical solution of one single peakon solution: (a) t = 1.0; (b) t = 2.0; (c) t = 3.0; 
(c) t = 4.0. 



Table 1: Comparison of and 7i errors for one-soliton propagation. 



At 


T 


Loo 




El 




MFE 0.001 


2.0 


6.9(- 


-4) 


8.1(- 


-4) 


0.001 


4.0 


1.4(- 


-3) 


1.6(- 


-3) 


RK4 0.01 


2.0 


1.7(- 


-12) 


4.5(- 


-13) 


0.01 


4.0 


5.4(- 


-12) 


8.4(- 


-13) 
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2r 

1.8- 
1.6- 
1.4- 
1.2- 
i 1 - 
0.8 
0.6 
0.4- 
0.2 



11 M 



(a) 



40 60 
X 



2 
1.8 
1.6 
1.4 
1.2 
< 1 
0.8 
0.6 
0.4 
0.2 
0, 



(b) 



40 60 80 

X 





2 
1.8 

1.6 
1.4 
1.2 
< 1 
0.8 
0.6 
0.4 
0.2 
0, 







1 i 


\ 
\ 
1 

V 



Figure 5: Numerical solution for collision of two nearly-peakon with pi ~ 198.9975, P2 = 
199.4995 and c = 200.0: (a) t = 0.0; (b) t = 10.0; (c) t = 15.0; (d) t = 20.0; (e) t = 30.0. 
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Table 2: Loo and Ii errors for two approximate peakon interaction by the self-adaptive moving 
mesh method 



At 


T 


Loo 




El 




MFE 0.001 
0.001 


5.0 
10.0 


2.2(- 
7.1(- 


-2) 

-2) 


5.5(- 
1.2(- 


-3) 
-2) 


RK4 0.01 
0.01 


5.0 
10.0 


2.0(- 
3.2(- 


-9) 
-9) 


1.5(- 
1.4(- 


-7) 
-5) 



1. Due to the integrability of the scheme and the self-adaptive feature of the 
non-uniform mesh, the Lgo-norm is small and the first conserved quantity is 
preserved extremely well even for a small number of grid points. 

2. The errors is mainly due to the time advancing methods. MFE is first order 
in time, so it produces relatively large L^o and Ei, roughly changing in pro- 
portional with time. RK4 is fourth-order in time, so up to T = 4.0, L^o and 
El are of the orders 10^^^ and 10^^^ for a grid number of = 101 and a 
time step At = 0.01. 

4.2. Propagation and interaction ofcuspon solutions 

The classical 4th-order Runge-Kutta method fails whenever the cuspon solu- 
tion is involved. It seems that a kind of instability occurs in this case, whose 
theoretical reason is still unclear. Therefore, only MFE is employed to conduct 
the numerical experiments hereafter. 

Example 3: One-cuspon propagation. The parameters taken for the one- 
cuspon solution are p = 10.98, c = 10.0. The number of grid is taken as 101 in 
an interval of width of 4 in the x-domain. Through the hodograph transformation, 
this corresponds to an interval of width 74.34 in the X-domain. Figure[6ta) shows 
the initial profile and the initial mesh. Figures [6tb)-(d) display the numerical 
solutions (solid line) and exact solutions (dotted line) at t = 2, 3, 4, together with 
the self-adjusted mesh. It can be seen that the non-uniform mesh is dense around 
the cuspon point, and moves to the left in accordance with the movement of the 
cuspon point. Table [3] exhibits the results of relative errors in Loo-norm and Ei. 

Example 4: Two-cuspon interaction. The parameters taken for the two- 
cuspon solutions are = 11.0, p2 = 10.5, c = 10.0. Figures |7ta)-(d) display 
the process of collision at several different times, along with the exact solution. 
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Table 3: Relative errors in Loo norm and the first conservative quantity for one-cuspon propagation 



At 


t 


Loo 




El 




0.005 
0.005 


2.0 
4.0 


3.3(- 
9.7(- 


-2) 
-2) 


4.7(- 
1.2(- 


-2) 
-1) 


0.001 
0.001 


2.0 
4.0 


1.1(- 
2.9(- 


-2) 
-2) 


1.2(- 
3.7(- 


-2) 
-2) 




Figure 6: Numerical solution of one single cuspon solution: (a) t = 0.0; (b) t = 2.0; (c) t = 3.0; 
(c) t = 4.0. 
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Meanwhile, the self-adaptive mesh is also shown in the graph. It can be seen that 
two cuspon solutions undertake elastic collision, regaining their shapes after the 
collision is complete. As mentioned in flS*], the two cuspon points are always 
present during the collision. The grid points are automatically adapted with the 
movement of the cuspons, and are always concentrated at the cuspon points. In 
compared with the exact solutions, we can comment that the numerical solutions 
are in a good agreement with exact solutions. As far as we know, what is shown 
here is the first numerical demonstration for the cuspon-cuspon interaction. 




10 20 30 40 50 60 10 20 30 40 50 60 

X X 



Figure 7: Numerical solution for the collision of two-cuspon solution with pi — 11.0, P2 — 10.5, 
c 10.0: (a) t = 13.0; (b) t = 14.8; (c) t = 16.6; (d) t = 25.0. 
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4.3. Soliton-cuspon interactions 

Here we show two examples for the soliton-cuspon interaction with c = 10.0. 
In FigjSl we plot the interaction process between a soliton of pi = 9.12 and a 
cuspon of p2 = 10.98 at several different times where the soliton and the cuspon 
have almost the same amplitude. It can be seen that when the collision starts 
(t = 12.0), another singularity point with infinite derivative (Wx) occurs. As 
collision goes on (t = 14.4, 14.6, 14.8), the soliton seems 'eats up' the cuspon, 
and the profile looks like a complete elevation. However, the cuspon point exists 
at all times, especially, at t = 14.6, the profile becomes one symmetrical hump 
with a cuspon point in the middle of the hump. 

In Figl9l we present another example of a collision between a soliton (pi = 
9.12) and a cuspon {p2 = 10.5) where the cuspon has a larger amplitude (2.0) 
than the soliton (1.0). Again, when the collision starts, another singularity point 
appears. As collision goes on, the soliton is gradually absorbed by the cuspon. 
At t = 10.3, the whole profile looks like a single cuspon when the soliton is 
completely absorbed. Later on, the soliton reappears from the right until t = 16, 
the soliton and cuspon recover their original shapes except for a phase shift when 
the collision is complete. 

4.4. Non-exact initial value problems 

Here, we show that the integrable scheme can also be applied for the initial 
value problem starting with non exact solutions. To the end, we choose an initial 
condition whose mesh size is determined by 

5k = 2ch{l-0.8sech{2kh-Wx/2)), (22) 

then, the initial profile can be calculated through the second equation of the semi- 
discretization, which is plotted in Fig|10](a). Figures [10] (b), (c) and (d) show the 
evolutions att = 10, 20, 30, respectively. Note that c = 10 in this computation. It 
can be seen that a soliton with large amplitude is firstly developed, and moving fast 
to the right. By t = 30, a second soliton with small amplitude is to be developed. 

Next, we increase the value of c to 90, which implies a very small dispersion 
term, corresponding to the dispersionless CH equation. The initial profile and the 
evolutions at t = 50, 150, 200 are shown in Fig{TT] It is seen that four nearly- 
peakons are developed from the initial profile at t = 50. Later on, an array of 
nearly-peakons of seven and eight are developed at t = 150, 200, respectively. 
This result is similar to the result for the KdV type equations with a small disper- 
sion, i.e. the peakon trains are generated. (For the KdV type equations, soliton 
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Figure 8: Numerical solution for cuspon-soliton collision with pi = 9.12, p2 = 10.98 and c = 
10.0: (a) t = 0.0; (b) t = 12.0; (c) t = 14.4; (d) t = 14.6; (e) t = 14.8; (f) t = 17.0; (g) t = 25.0. 
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Figure 9: Numerical solution for cuspon-soliton collision with pi = 9.12, p2 = 10.5andc = 10.0: 
(a) t = 0.0; (b) t = 9.0; (c) t = 10.0; (d) t = 10.3; (e) t = 10.6; (f) t = 11.5; (g) t = 16.0. 
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Figure 10: Numerical solution starting from an initial condition (19) with c = 10: (a) t = 0.0; (b) 
t = 10.0; (c) t = 20.0; (d) t = 30.0. 
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trains are generated. For example, see [|44l . l45h for numerical simulations and 11461] 
for a theoretical analysis for the KdV equation.) A theoretical analysis for the dis- 
persionless CH equation to explain the above intriguing numerical result is called 
for. 
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Figure 11: Numerical solution starting from an initial condition (19) with c — 90: (a) t ~ 0.0; (b) 
t = 50.0; (c) t = 150.0; (d) i ^ 200.0. 



5. Concluding Remarks 

In the present paper, we have proposed a self-adaptive moving mesh method 
for the CH equation, which based on an integrable semi-discretization of the CH 
equation. It has the properties: (1) it is integrable in the sense that the scheme 
itself admits the A^-soliton solution approaching to the iV- soli ton solution of the 
CH equation in the limit of mesh size going to zero; (2) the mesh is non-uniform 
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and is automatically adjusted so that it is concentrated in the region where the so- 
lution changed sharply, for example, the cuspon point; (3) once the non-uniform 
mesh is evolved, the solution is determined from the evolved mesh by solving a 
tridiagonal linear system. Therefore, either from the accuracy or from the com- 
putation cost, the proposed method is expected to be superior than other existing 
numerical methods of the CH equation. This is indeed true. The numerical results 
in this paper indicate that a very good accuracy is obtained. 

Two time advancing methods, the modified forward Euler method and the 
classical 4th-order Runge-Kutta method, are used to solve the evolution of non- 
uniform mesh. The Runge-Kutta method gains much better accuracy than the 
modified forward Euler method. However, it fails for the computations of cuspons. 
Using the self-adaptive moving mesh method for the CH equation, we have ob- 
tained interesting numerical computation results starting with non-exact solutions. 
When K is very small, the peakon train is generated from the non-exact initial con- 
dition. 

As further topics, it is interesting to construct integrable discretizations, or, 
the self-adaptive moving mesh methods for a class of integrable nonlinear wave 
equations possessing soliton solutions with singularities such as peakon, cuspon 
or loop solutions. For example, such equations include the short pulse equation 
which was derived as a model for the propagation of ultra-short optical pulses in 



nonlinear media [14711 . 

uxT = u + -{u^)xx, (23) 

D 



and the Degasperis-Procesi (DP) equation 114811 

ut + Sk^ux - utxx + 4mmx = ^uxUxx + uuxxx- (24) 

It is worth pointing out that the authors have constructed semi- and full-discretization 
for the short pulse equation, which is another example of the self-adaptive mov- 
ing mesh method (a9], in which we have succeeded in computing the one- and 
two-loop soliton propagations and interactions. 
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